Canonical-basis solution of the Hartree-Fock-Bogoliubov equation on 
three-dimensional Cartesian mesh 
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A method is presented to obtain the canonical-form solutions of the HFB equation for atomic 
nuclei with zero-range interactions like the Skyrme force. It is appropriate to describe pairing 
correlations in the continuum in coordinate-space representations. An improved gradient method is 
used for faster convergences under constraint of orthogonality between orbitals. To prevent high- 
lying orbitals to shrink into a spatial point, a repulsive momentum dependent force is introduced, 
which turns out to unveil the nature of high-lying canonical-basis orbitals. The asymptotic properties 
at large radius and the relation with quasiparticle states are discussed for the obtained canonical 
basis. 

PACS numbers: 21.10.Pc, 21.30.Fe, 21.60.Jz, 27.30. -(-t 



I. INTRODUCTION 

Pairing correlations play an essential role in the deter- 
mination of the ground-state structure of the vast ma- 
jority of atomic nuclei. Among its treatments, an easy 
but still enough accurate one is to consider only single- 
particle states near the Fermi level while the effects of the 
other states are assumed to be absorbed in the strength 
of an effective pairing interaction. In such a treatment, 
our experiences in Hartree-Fock(HF)-|-BCS calculations 
suggest that at least half of a major shell above the Fermi 
level must be considered for meaningful estimations. As 
a consequence, one has to explicitly consider positive- 
energy states if the Fermi level of neutrons is higher than 
the negative of half of the major shell spacing (~ ^A^^/"^ 
MeV). This condition applies to about half of the 10"' nu- 
clides which exist between proton and neutron drip lines 
in the nuclear chart, not only to nuclei near the neu- 
tron drip line or outside the s-process path. Thus there 
are plenty of necessities for theoretical frameworks which 
can describe pairing correlations involving the continuum 
states. 

If the pairing correlation significantly involves the con- 
tinuum states, the HF-I-BCS approximation becomes in- 
adequate because the occupation of unbound HF orbitals 
leads to the unwanted dislocalization of the nucleon den- 
sity. This is a serious problem in coordinate-space treat- 
ments, which is more favorable to describe loosely bound 
systems like drip-line nuclei than expansions in harmonic- 
oscillator eigenstates (except the transformed oscillator 
basis of Ref. fl). The reason why a dislocalized solu- 
tion becomes the variational minimum is that, in order 
to separate the variational equations into HF and BCS, 
one has to neglect the effects that the matrix elements of 
pair-scattering processes are affected by the changes in 
the wave functions of the orbitals involved 0. To fully 
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take into account these effects leads to the Hartree-Fock- 
Bogoliubov (HFB) equation, with which the density is 
localized whenever the Fermi levels are negative 0, Q| . 

The HFB method in the coordinate space was first 
formulated using the quasiparticle states and solved for 
spherically symmetric states in Ref. . Although spher- 
ical solutions can be obtained easily with present com- 
puters (for zero-range forces), deformed solutions are still 
difficult to obtain because there are quite a large num- 
ber of quasiparticle states even for a moderate size of the 
normalization box (i.e., the cavity to confine the nucleons 
to discretize the positive-energy single-particle states). 
An orthodox approach to face this difficulty is the two- 
basis method IE S S 13 ' which the quasiparticle 
states are expanded in bound and unbound HF orbitals. 
This method requires heavy numerical calculations be- 
cause there are a pile of unbound HF orbitals below an 
energy cut-off of even only a few MeV. 

An alternative approach is the canonical-basis HFB 
method. According to the Bloch-Messiah theorem [lol| . 
every HFB solution obtained as the vacuum of a set of 
Bogoliubov quasiparticles has an equivalent expression in 
terms of a BCS-type wave function. The single-particle 
states appearing in this expression are called the HFB 
canonical basis. In the canonical-basis HFB method, one 
obtains the solution in the canonical form without using 
the quasiparticle states. This method appeared originally 
in Ref. llj to obtain spherical solutions. However, there 
are no sever difficulties in obtaining HFB solutions for 

)ds. Appfi catic 

a three-dimensional Cartesian mesh representation 15]. 
Incidentally, a different line of application is also found 
in literature [T^ . 

Let us explain the advantage of the canonical-basis 
method over the two-basis method concerning the treat- 
ment of the pairing in the continuum using Fig. ^ On 
both sides of the figure, the ordinate represents the ex- 
pectation value of the mean-field (HF) energy, while the 
abscissa stands for the radius r from the center of the 
nucleus, ec and ep mean a cut-off energy and the Fermi 



spherical nuclei using any other methods. Appli cations to 
deformed nuclei have been done by us |3. ll2l[l3llT^ u sing 
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FIG. 1: Comparison between the HF orbitals and the HFB 
canonical-basis orbitals. See text for explanations. 



level, respectively. On the left-hand side, the wavy hori- 
zontal lines stand for energy levels and their spatial ex- 
tent for the HF potential denoted by a solid curve. The 
wave functions of positive energy states extend to the wall 
of the box and their level density is much larger than that 
of negative energy states. In the two-basis HFB method, 
one has to mix these positive-energy orbitals to construct 
localized canonical-basis orbitals, which is a numerically 
demanding task. On the right-hand side, the wavy lines 
represent the HFB canonical-basis orbitals. Unlike HF 
orbitals, they are spatially localized for both negative and 
positive energies. Because of this localization, the level 
density is much smaller than the unbound HF orbitals. 
Therefore, one needs much less orbitals. More specifi- 
cally, the number of necessary single-particle states to 
describe the HFB ground state of a nucleus is propor- 
tional to the volume of the nucleus in the canonical-basis 
method while it is related to the volume of the normal- 
ization box in the other methods. Incidentally, with the 
dash curve, we suggest the existence of some potential 
which binds the high- lying canonical-basis orbitals. The 
identity of this potential is unveiled in Sec. IVIBI 

In this paper we formulate the canonical-basis HFB 
method on a cubic mesh and develop an efficient gradient 
method to obtain its solutions. In order to decrease un- 
physical influences of high-momentum components due 
to zero-range interactions, we introduce a momentum de- 
pendent term to the pairing interaction and show how it 
suppresses a problematic behavior of wave functions pe- 
culiar to the canonical-basis HFB method. As a byprod- 
uct, we show a clear way to understand the nature of 
high- lying canonical-basis orbitals. 

The contents of this paper is as follows. In sections 
m ~ E we give the formulation and general considera- 
tions. Discussions using results of numerical calculations 
are collected in section IVII Conclusions are given in sec- 
tion In Appendix 1X1 we examine the errors due to 
the mesh representation and the box boundary condition. 
Appendix IbI gives a model analysis of the problem pecu- 
liar to the method, i.e., a phenomenon that canonical- 
basis orbitals collapse to a spatial point once their occu- 
pation probabilities fall below some critical value when 
pure delta- function forces are used. 

Some of the contents of this paper have already been 
discussed in our earlier reports 0, 0, 0, 0| . For ex- 
ample, Ref. includes the first report on the canonical- 
basis HFB method on a cubic mesh. Refs. ^3 



ported on the discovery of the point-collapse problem 
and the tests of the state-dependent cut-off factors and 
pairing-density-dependent forces as the remedies. In this 
paper we have rewritten those contents to include the 
momentum-dependent interaction and, by including the 
rewritten contents, organized the formalism for consis- 
tency and selfcontainedness. All the numerical calcula- 
tions shown in this paper have been performed using the 
renewed form of the interaction. 

Incidentally, some important parts of those reports are 
not repeated in this paper. The examples are an analysis 
of the spherical two-basis method concerning the preci- 
sion of density tails @ and numerical demonstrations of 
the point collapses of high- lying orbitals Q . 



II. HFB IN THE CANONICAL 
REPRESENTATION 

To begin with, let us formulate HF and HFB methods 
in the coordinate-space representation in order to eluci- 
date the difficulty of the quasiparticle-based HFB method 
compared with the HF method and to propose its possi- 
ble solution in terms of the canonical-basis HFB method. 
For the sake of simplicity, we consider only one kind of 
nucleons in this section. The number of that kind of nu- 
cleons are designated by A''. The ^-component of the spin 
of a nucleon is represented by s (= ±i). 

In HF, one should minimize (\l/|iJ|\E') among all the 
A^-body single Slater-determinant states, 

N 
i=l 

°I = E / d'ri:,{r,s)a\r,s), (2) 

s •' 

by varying {4'i}i=i,---,N under orthonormality conditions 
{ipi\ipj) = 5ij. The operator a| creates a nucleon with 
a wave function ipi{r,s). The ket state |0) stands for 
the state in which no nucleons exist. The distribution 
function of the number density of the nucleons is related 
to the single-particle wave functions as 

N 
i—l s 

There is arbitrariness in the selection of {ipi} as for uni- 
tary transformations among them, tpi are called HF or- 
bitals when they diagonalize the HF Hamiltonian. 
In HFB, the solution takes the following form. 
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where bi is the annihilation operator of a positive-energy 
Bogohubov quasiparticle with two types of ampUtudes 
(j)i{r,s) and ipi{r,s) for particle-Hke and hole-hke parts 
of the excitation. A^b is the number of the basis states 
of the representation used and is by far larger than N. 
One should vary (^i}i=i,...,AfB under an orthonormal- 
ity condition 

(6) 

and a constraint on the expectation value of the number 
of nucleons, 



p{r)d^r = N, 



(7) 



where the nucleon density for state Q is expressed as 



The essential difference between HF and HFB lies in 
the number of necessary single-particle wave functions. 
It is N in the former and 2A'b in the latter. Obviously, 
the latter is much larger. 

Owing to the Bloch-Messiah theorem , the state (@J) 
can be transformed (as for the ordinary ground states of 
even-even nuclei) into the following form, 



I*) - n( 



Ui + Vi a] a\ 



|0), 



E 



(9) 



(10) 



where a\ and a\ create a nucleon with wave functions 
tpi{r,s) and ^i{r,s), respectively, which are called the 
canonical basis of the HFB vacuum j^*). They may 
be called the natural orbitals^ instead, because they 
are the eigenstates of the one-body density matrix. In 
this paper, we often call them canonical-basis orbitals or 
canonical orbitals. It is required N^j = ^N-q for the exact 
equivalence between Eqs. I^J and Q to hold in general. 

In this paper we assume that I^E*) is invariant under the 
time reversal operation. In this case, V'i and "04 are the 
time-reversal state to each other and only one of them has 
to be considered explicitly in the variational procedure. 
The nucleon density can be expressed as 



(11) 



in the time-reversal invariant case. 

In order to obtain solutions expressed in the canonical 
form without using the information of the quasiparticle 
states, one should vary {V^i, u^, Wi}i=i^..._jv„ under three 
kinds of constraints, i.c, the orthonormality condition, 



(12) 



a condition on the expectation value of the number of 
nucleons, 



2E^ 



(13) 



and the normalization of the u-v factors, 



vf = 1. 



Reinhard et al. wrote that the advantage of the repre- 
sentation @ over ^ is that one has to consider only a 
sing le set of wave functions {ipi}, not a double set ipi} 
|ll|. In our opinion, however, one may expect much 
greater benefit from the canonical-basis representation. 
Namely, i may be truncated as i < N„ = 0{N) ^ 
to a very good approximation. 

The reason is the localization of the density. The den- 
sity distribution of HFB solutions can be shown to behave 
asymptotically for large r as 



p(r) 



K > 



(14) 



as far as the Fermi level ep is negative Q , where m is the 
nucleon mass. Consequently, ipi appearing on the right- 
hand side of Eq. Hll|l must be a localized function as p{r) 
on the left-hand side. On the other hand, the orthogonal- 
ity relation (|12|) restricts the number of wave functions 
which can exist in the vicinity of the nucleus. Therefore 
the number of canonical orbitals cannot be very large. In 
mesh representations, A^b is proportional to the volume 
of the box while N^j is proportional to the volume of the 
nucleus. The latter is 10~^ — 10^^ times as small as the 
former in typical calculations. 

Situation is quite different in the quasiparticle formal- 
ism: The number of quasiparticle states is proportional 
to the volume of the box. This is because the localization 
of the density demands only that of (pi through Eq. (jH)) 
while (pi does not have to be localized in general. The 
orthogonality condition © involves both ipi and 0i . This 
enables many quasiparticle states having similar tpi to be 
orthogonal to each other by differing their (pi. 



III. MEAN FIELDS FOR ZERO-RANGE 
INTERACTIONS 

In this paper, we employ the Skyrme interaction |3.IT^. 
which is a density- and momentum-dependent zero-range 
interaction: 

+ t2{l+X2Pa)k-dk+lt3{l+X3Pa)p"S, (15) 

where Pa is the spin exchange operator, S — S {ri — 
with ri and r2 the position vectors of the interacting two 
nucleons, p is the nucleon density at ri (= r2), and 



(16) 
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Tik is the relative momentum between the two nucleons. 
Since it is hermitian under the box (vanishing and peri- 
odic) boundary conditions, we have not specified in which 
way (to bra or ket states) it operates in Eq. IjlSfl unlike in, 
e.g., Ref. Zero-rangeness makes the mean- field poten- 
tials local, which is an essential advantage for coordinate- 
space solutions. 

Among the terms of the Skyrme force, those with 
strength parameters to, ti, and ^3 act only on s-wave 
relative motions. The term with ti serves to take into 
account the effects of the finite-rangeness in terms of 
the lowest order momentum dependence, while the term 
with expresses a density dependence. The term with 
strength t2 acts on relative p-wave states through a 
different form of momentum dependence. Note that 
the isospin channel (T^O or 1) is uniquely determined 
through the requirement of the two-body antisymmetry. 
Owing to the p-wave as well as the s-wave interaction 
terms, the Skyrme force can act on all the four kinds of 
spin-isospin channels. 

The complete Skyrme interaction includes a spin-orbit 
term iW{(Ti + ar^) ■ k x Sk, which we have excluded from 
Eq. H15|) . This elimination decreases the size of the nu- 
merical calculations by a factor of four (two from spin 
up and down components, another two from real and 
imaginary parts of the spatial wave function). Due to 
this size reduction, we can perform very severe tests of 
the canonical-basis HFB method by, e.g., taking into ac- 
count unusually many canonical orbitals or expressing 
wave functions on a very large mesh. 

As the parameters of the force, we choose the SIII set 
[r^ except for the spin-orbit term, i.e., to = —1128.75 
MeV fm^, <3 =14000 MeV fm^+S", a = 1, and = 0. As 
for the parameters appearing in the momentum depen- 
dent terms, the following two combinations are sufficient 
to treat N — Z even-even nuclei. 



/l = j^(3ii+5i2+4i22:2), 
/2 = ^(-9tl + 5i2+4t2.T2), 



(17) 
(18) 



whose values are /i — 44.375 MeV fm^ and /2 = —62.969 
MeV fm^. 

As is usually done, we employ different parameters be- 
tween the mean- field and the pairing interactions. We 
use the following interaction for the pairing channel. 



P 2 



l-^-f^l U 



(19) 



where Vp is the strength and ^(1 — P^) is the projector 
into spin-singlet states. Terms in the braces represent 
dependences on particle and pairing densities. The in- 
teraction becomes repulsive where p > pc or p > pc- For 
the particle-density dependence, we use pc = 0.32 fm~'^ to 
vanish the volume-changing effect (This choice, twice 
of the matter density, is also recommended for a different 



reason |20j.) The dependence on the pairing density p (to 
be defined by Eq. ij^HI)) has been introduced in Ref. [l3| 
to prevent unphysical behaviors of high-lying canonical 
orbitals discussed in Appendix IbI This term is squared 
because p can become negative in principle. Except in 
Appendix \n[ we use pc = 00. The momentum depen- 
dent term acts on s-wave relative motions and quenches 
interactions between nucleons in high-momentum states. 
The interaction vanishes at relative momentum kc at low- 
density points. If Pc = kc = 00, this pairing force is re- 
duced to that introduced in [13. The typical values of 
the parameters are Vp = —1000 MeV fm"^ and kc=2 im~^ 
but different values are also used. 

The S^O and T=l part of the Skyrme force of Eq. ^ 
may be rewritten as Eq. (|19l) without the pairing density 
dependent term. Using different parameterizations be- 
tween the two forces has an advantage of avoiding con- 
fusions. 

When one considers both protons and neutrons, the 
state of the nucleus is usually assumed to be a product 
of two BCS-typc states of Eq. ©, one for the protons 
and the other for the neutrons: 



1*) = 



q=p 1=1 



"qi "qi 



|0), 



(20) 



where the index q distinguishes between protons (p) and 
neutrons (n). a^^ creates a proton having a wave function 

'0pi(r, s) while aj^,- creates a neutron with a wave function 
ipnii'i'Ts). This product form matches our choice of the 
pairing interaction of Eq. (|19|l which acts only on T = 1 
pairs. 

In this paper we treat only N=Z nuclei without 
Coulomb interaction for the sake of simplicity. In this 
case, the wave functions are the same between protons 
and neutrons. Moreover, because the potentials are in- 
dependent of the spin, the wave functions 'ipi{r,s) can 
be factorized into a product of a spin wave function and 
a real-number function of the position, which we write 
i/jiir) hereafter. More specifically, we assume ippi{r,s) 
= ipni{r,s) = 'ipi{r)S^ i and i^pi{r,s) = i}niir,s) = 

The isoscalar pairing may take over the isovcctor pair- 
ing in ~ Z nuclei |21j. Nevertheless, we assume 
Eq. (|20|l because the aim of this paper is to examine 
a method whose most important applications are to the 
nuclei near the neutron drip line. 

Our Hamiltonian H consists of a kinetic energy term 
and two-body interaction terms. It should be understood 
that both are expressed in the second quantized form be- 
cause HFB states do not have a fixed number of particles. 
The kinetic energy term is —h^W"^ /2m in the one-body 
expression, for which we assume the mass of a nucleon 
m to be the average of the proton and the neutron bare 
masses (nip and mn) with an approximate correction fac- 
tor for the center-of-mass motion of a nucleus of mass 
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number A: 



(21) 



For the interaction terms, we use the interaction (|15() for 
the mean-field type contractions and the interaction H19|) 
for the pairing type contractions in evaluating the ma- 
trix elements of the two-body interactions using Wick's 
theorem. 

The total energy for the state H20() can be expressed in 
terms of a single space integral as, 



(22) 



where the integrand is given by 



top + fipr + /2PV V 




16 



2+a 



.(23) 



On the right-hand side of the above equation, functions 
of position r are 



p{r) 
r(r) 
p{r) 
f(r) 



i=l 

i=l 

4^u,i;,|VV'»(r-)|^ 



(24) 
(25) 
(26) 
(27) 



i=l 



which are called the (particle) density, the kinetic- 
energy density, the pairing density, and the pairing 
kinetic-energy density, respectively. 

There is a set of single-particle operators (Ti-i, h, and 
h) which are useful to obtain mean-field solutions. They 
can be obtained by taking the functional derivative of 
^'tot with respect to a canonical wave function ipi. This 
turns out to be equivalent to operating a state-dependent 
single-particle Hamiltonian Tii to the wave function : 



SEu 



Hi 



(28) 



Here, we have taken the familiar point of view that ipi and 
ip* can be formally regarded as independent variables, 
although ipi is actually a real function. 

One can express Tii, in turn, as a linear combination 
of state-independent single-particle Hamiltonians as. 



Til = v,h + UiVih, 



(29) 



where h and h are called the mean-field and pairing 
Hamiltonians. They are given by 



with 
B 
V 
B 



h = -V-BV + V, 
h = -V-SV + F, 



Itop + ^i3p'+" + /it + 2/2 V^p 



8fc2 



V 



Pc 



Pc 



2fc? + 



(30) 
(31) 



(32) 

,(33) 
(34) 

.(35) 



We call V and V the mean-field and the pairing poten- 
tials, respectively. 

When there are no pairing correlations, the mean- 
field Hamiltonian h is identical to the HF single-particle 
Hamiltonian. Therefore, its expectation value. 



{ipi\h\i>i 



(36) 



may be called the single-particle energy of the ith canon- 
ical orbital. The negative of the expectation value of the 
pairing Hamiltonian, 



A. 



(37) 



has a meaning of the pairing gap of the orbital. The 
pairing gap should be positive or zero. Indeed, could 
take negative values only if the wave function would con- 
sist of very high-momentum components . A related 
discussion is given in Sec. IVIBI 

The quasiparticle states of Eq. (jsj are the eigenstates 
of their own Hamiltonian, i.e., so-called HFB super ma- 
trix composed of h and h: 



h — ep 
h 



qp 



(38) 



where ep is the Fermi level. In small subspaces like sev- 
eral major shells, the easiest way to obtain an HFB so- 
lution is to solve Eq. (|38|l as an eigenvalue problem of a 
hermitian matrix. In large spaces like mesh representa- 
tions, however, it is easier to obtain the canonical orbitals 
directly, rather than via quasiparticle states. Such a di- 
rect method is explained in the next section. 

Incidentally, if one considers seriously that ipi is a real 
function, one may notice that the functional derivative 
should be doubled since 



6Et, 



SEtot Sip. 



bra 



6Etot Sip. 

kct ' 



kot 



2 • AH^^P^, (39) 



as V-'i"'* = "01 and ipf^^ = Eq. (|39|l seems to be a 
more appropriate expression of what our numerical cal- 
culation program actually does. However, the factor 2 at 
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the beginning of the right-hand side as well as the factor 
4 which reflects the spin-isospin degeneracy do not mat- 
ter because they can be conveniently absorbed into the 
step-size parameter of the gradient method (see Sec. 01. 



IV. CUT-OFF SCHEMES FOR THE PAIRING 
INTERACTION 



In this section we discuss on the cut-off schemes 
for zero-range pairing interactions in relation to the 
canonical-basis HFB method. 

Delta function forces without a cut-off lead to a diver- 
gence of the strength of the pairing correlation. In order 
to prevent the divergence, one usually takes only those 
quasiparticles whose excitation energy Eqp is lower than 
some cut-off energy Ec in constructing the ground state 
0. Namely Eq. igj) is modified to 



I*) = 



qp 



6.10), 



(40) 



where 9 is the Heaviside (step) function. 

In this paper we do not discuss how one should change 
the strength as a function of Ec for renormalization. In- 
stead, we regard Ec as one of the constants to define the 
force. 

In order to implement the cut-off, one should pay at- 
tention to the following two peculiarities of the canonical- 
basis HFB method, i) The introduction of Ec causes a 
serious inconvenience to this method, i.e., the absence 
of the quasiparticle Hamiltonian. For this reason, it is 
preferable to use the number of canonical orbitals iV^ in- 
stead of Ec- ii) For delta- function pairing interactions, 
substitution of Ec with leads to a physically meaning- 
less situation, in which canonical orbitals are collapsed to 
a spatial point (a point collapse phenomenon, hereafter). 
To overcome this difficulty, we modify the pairing inter- 
action by adding a repulsive momentum-dependent term. 
In the rest of this section we discuss these two problems. 

i) Let us explain how the introduction of the cut-off 
causes the absence of the quasiparticle Hamiltonian. An 
analogy with the BCS approximation indicates a natural 
way to introduce an energy cut-off to the canonical-basis 
method. Namely, the smooth cut-off method [Tsll modifies 
the seniority interaction as 




/ \j>0 / 



(41) 



where fi = f{ei) is a function of single-particle energy 
Ei and takes on ~ 1 well below a chosen cut-off energy 
and smoothly becomes zero above it. The smoothness is 
indispensable to the stability of the solution. In analogy, 
we may modify the pairing density as 



= 4^ UiVi lipi 



4^ /i UiVi IV-i 



(42) 



with 

fi = exp 



ip*ir)V^^P,{r)d^r. (43) 



In Eq. (|43|) . we assume a dependence on the kinetic 
energy (cx kf), not on the mean-field energy Ci as in 
Eq. H41|) , to avoid a highly complicated expression of the 
gradient for the latter case. (In BCS, such complications 
are simply neglected.) The gradient is given for kc = oo 
by 



1 -^^tot 

4 S^P* 



vfh + u,v,h{f{kf) + f{kf)V^} A- (44) 



This method works well to prevent the point collapses in 
numerical tests using /i = 1.2 fm taken from the Gogny 
force [23. However, this cut-off scheme has an disadvan- 
tage that the HFB super matrix in Eq. H38|) cannot be 
defined, because in Eq. H44|l the pairing Hamiltonian (the 
part proportional to UiVi in square brackets) is dependent 
on the canonical-orbital index i. This means that there is 
no common pairing Hamiltonian to all the quasiparticles. 
Consequently one does not have well-defined Bogoliubov 
quasiparticle states and cannot utilize them to refine the 
HFB solution, which is a rather serious drawback as dis- 
cussed in the next section. If one can use in place 
of Ec, however, one can recover a unique quasiparticle 
Hamiltonian. 

Incidentally, a similar kind of ambiguity arises from 
the cut-off scheme in terms of the quasiparticle energies 
according to Eq. H40|) . In this case, the canonical orbitals 
do not have well-defined Hamiltonians. It may not mat- 
ter, however, because one does not need Hamiltonians to 
define canonical orbitals when they can be obtained by 
transforming the quasiparticle states. 

ii) Let us discuss the phenomenon of the point collapse. 
In quasiparticle (or two-basis) HFB method, specifying 
the number of quasiparticle (HF single-particle) states to 
be considered has the same consequences as imposing an 
energy cut-off because the states to be chosen are anyway 
the eigenstates of the lowest energies of the quasiparticle 
(HF) Hamiltonian. 

On the other hand, the canonical-basis HFB method 
is very peculiar concerning this point. The canonical or- 
bitals are not the eigenstates of a single operator but each 
of them has its own Hamiltonian Tii as a function of the 
occupation probability given by Eq. (|29() . Hi becomes h 
and h in the two extreme situations vf=l and 0, respec- 
tively. It is only when one considers a small number of 
orbitals that all of them are the lowest-energy approxi- 
mate eigenstates of h. Otherwise, some of the orbitals 
may become the eigenstates of h rather than h: We show 
in Appendix IbI that, if the occupation probability of an 
orbital becomes sufficiently small, it spontaneously de- 
creases the probability further toward zero (to decrease 
the total energy of the nucleus) and change its Hamilto- 
nian to h. This occurs when the pairing force is a delta 
function (with or without particle-density dependence). 
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for which h is a local potential and thus its eigenstates 
are delta functions in the spatial coordinate. This means 
that their (h) is infinitely high. On the other hand, the 
rest of the orbitals are left in low-energy subspacc of h. 
Therefore the restriction on Ny^ results in utterly different 
solutions from those obtained by using a cut-off energy 
explicitly. 

This peculiarity of the canonical-basis method is essen- 
tial to the implementation of the method. In the quasi- 
particle and the two-basis methods, the unphysical diver- 
gence of the pairing correlation strength occur only in the 
limit of oo (supposing 7V„ is used instead of Ec). 

In the canonical-basis method, in addition to this diver- 
gence, a different type of unphysical behavior (the point 
collapse) can also happen for fixed at finite numbers. 
One has to take special care of the latter problem if one 
tries to implement the canonical-basis method. 

Now, if one can suppress the point collapse, one is al- 
lowed to use iVw as the control parameter of the cut-off 
energy. Appendix ^ shows that this can be achieved, 
e.g., by introducing a sufficiently strong repulsive mo- 
mentum dependence to the delta-function pairing inter- 
action. With this term applied to finite (small or large) 
numbers of canonical orbitals, unphysical situations in- 
cluding the point collapse never happen. (If it is applied 
to infinite number of orbitals, not the point collapse but 
the divergence occurs 22J.) In this paper we employ this 
momentum-dependent force. 

In addition to the role to enable a cut-off in terms of 
iVw, equally important is the byproduct that it also clar- 
ifies the nature of high-lying canonical orbitals by pro- 
viding a kinetic energy term to the pairing Hamiltonian 
h. Further explanations using a numerical example are 
given m Sec. Em 

Incidentally, Appendix IbI also gives a discussion on an 
alternative way in terms of a pairing-density dependent 
force to enable cut-off in terms of Ny^ . 

Before ending this section, we mention the prospect 
about the renormalization of the pairing force strength. 
Considering the recent quantitative success of the regu- 
larization method of delta-function forces based on the 
Thomas- Fermi approximation |24,^5j 26j,(23|, we think it 
is worth trying in the canonical-basis method in future. 

However, it will require a special treatment. The 
Thomas-Fermi approximation makes the interaction 
strength a decreasing function of the particle density. 
One might expect this modification of the interaction 
would hinder the point collapse. Actually, dependences 
on particle density cannot prevent it because the col- 
lapse makes only the pairing density divergent but leaves 
the particle density unchanged (because v cx w^/^ and 
|V'w(0)P c!c as shown in Appendix Ib)| . Therefore one 
has to keep the momentum-dependent term. This will 
lead to a modification of the regularization procedure. 



V. GRADIENT METHOD FOR 
CANONICAL-BASIS HFB 

We explain in this section the details of a procedure to 
obtain the canonical-basis solution of the HFB equation 
directly, not by way of quasiparticle states. What should 
be done is to minimize iJtot of Eq. (|22|l under constraints 
of Eqs. H12|l and 1)13(1 . Equivalently, one may introduce a 
Routhian, 

R - Stot-4eFX]wf 

i=l 

N„ N„ 

- 4EE^^H(V'.|^j->-'5y}, (45) 

and minimize it without constraints, ep is probably the 
most familiar Lagrange multiplier, whose physical mean- 
ing is the Fermi level. In the definition (|45|) . Lagrange 
multipliers A^^- obeying hermiticity, 

A,, = A*„ (46) 

are introduced instead of ^Nvj{Nvj + 1) independent mul- 
tipliers. The hermiticity ensures the equality between 
the number of constraints and the number of indepen- 
dent multipliers. Moreover, it makes R real so that two 
conditions, SR/6ipi = and SR/Sip* = 0, become equiv- 
alent and thus one has to consider only one of them. 
Note that 6ij is subtracted from (ipilipj), in contrast to 
Ref. • This subtraction is necessary to treat A^ not as 
constants like ep but as functionals of the wave functions. 

The stationary condition of R results in two kinds of 
equations. One is dR/dvi — 0, which concerns the occu- 
pation amplitudes Vi and is fulfilled by [ll[ 

Among the double sign, the minus sign corresponds to 
the minimum. The relative sign between Vi and Ui = 
rk^/X — vl should be plus, i.e., UiVi > (unless < 0). 
The stationary condition for a wave function ^Ji{r) leads 
to the following equation: 

I SR ^ , , , , 

- EE^{(^.l^'=)-^.'=}-0' (48) 
j=i k=i 

Owing to the state dependence of Hi, the orthogonal- 
ity condition becomes nontrivial to fulfill. In HF, the 
orthogonality is automatically satisfied because are 
eigenstates of the same hermite operator h, which are 
orthogonal to one another. The orthogonalization proce- 
dure is needed only because of the instability for decaying 
into Pauli-forbidden configurations. On the other hand, 
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in the canonical-basis HFB method, the orthogonahza- 
tion is indispensable and the explicit functional form of 
Xij is the most important secret of the method. 



Reinhard et al. have proposed pj| 



(49) 



Let us give our reasoning on how the above form can 
be deduced. This consideration plays the crucial role in 
order to modify the form for the sake of a faster conver- 
gence. The requirement that Eq. (|48|l must hold at the 
solution (where (V'ilV'j) = 



means, 



"Hiltpi) = ^ AifclV'fc) (at the solution). (50) 



k=l 



By taking the overlap with {ipjl, one obtains, 
Xij = {ipj\T-Ci\ipi) (at the solution). 



(51) 



Eqs. H49I) and (|51l) are equivalent at the solution because 
Xij is defined to be hermite by Eq. H4t)|) . However, one 
should employ Eq. (|49|l rather than Eq. (|51() because the 
former, but not the latter, is hermite at points other than 
the solution. 

One can utilize the gradient method to obtain the HFB 
solutions in the canonical-basis formalism. Let us de- 
scribe this method briefly: Small variations in ipi and ijj* 
change the Routhian as, 



Sip, Sip* 



(52) 



which indicates the steepest descent direction to be, 
5R _ , 5R 



S'tpi oc 



Si'* 



Hi « 



Sipi 



(53) 



One takes a small distance movement toward this direc- 
tion, i.e., tpi ^ 4'i + S'ipi with 



S-ipi 



E 



(54) 



By repeating the evaluation of {Sipi} and the movement, 
one eventually reaches the minimum of R. Incidentally, 
in the braces on the right-hand side, the term including 
SXjk/Sip* in the second member of Eqs. (|48|l is dropped 
because we know empirically that the orthogonality is 
stable without this term. 

In Eq. H54() . At is a parameter introduced to control 
the size of a movement. One should regard that the prcf- 
actors on the right-hand side of Eq. (|39|) are included in 
this parameter. One may call At the imaginary-time- 
step size because the first order approximation in At of 
the imaginary-time evolution method becomes essen- 
tially the same procedure as the gradient method : 



Sip I 



e-^'^^i/j,-iP, 



Atn,ip, + O { At"^) 



(55) 



In order to justify the negli gen ce of the higher-order 
terms in Eq. H55|) . it must hold|2q| 



/cvo < 2 for At 



(56) 



where T^ax is the maximum kinetic energy (used as 
the upper bound of the single-particle energy). For the 
Cartesian mesh representation with a mesh spacing a. 



(57) 



where B{p) is given by Eq. %\2\ . By choosing p between 
p — (free space) and p—Q.lQ fm^^ (nuclear matter den- 
sity), one obtains Tmax = 1254 MeV for the SIII force 
with a = 0.8 fm~^. For the calculations shown in this 
paper, we have used /evo = l-O, i-e.. At = 5 x 10~^^ sec. 
Kinetic energy also arises from the pairing Hamiltonian 
h. However, an estimation by replacing B{p) of Eq. lfF7|l 
with B{p) of Eq. (|34|) results in an order of magnitude 
smaller value. 

What is presented above may be called a naive version 
of the gradient method. The convergence to the HFB so- 
lution turns out to be very slow with this naive method. 
A numerical example will be given in Sec. I VI ^ The ori- 
gin of this slowness can be traced back to the factor mvi 
in TLi, which is just a linear combination of h and h with 
coefficients w| and UiVi. The effects of Tii can be very 
weak for canonical orbitals having small Vi, which leads 
to the smallness of their changes in a gradient-method 
step. 

Now, steepest-descent paths depend on the choice of 
the variables. For example, Eq. (|54|l is obtained when one 
uses ipi as the variables. If one uses norm-resized wave 
functions Xi — V'i/v^: where ai is a scaling factor, a 
gradient-method step becomes Sxi oc —SR/Sx*, which is 
equivalent to 



Sipi 



At a., 



SR 



-ai At TLiipi - ^ Xijipj 



(58) 

Thus the path is changed. By choosing ai to cancel 
the small factors in TLi, one can accelerate the otherwise 
slow convergence, which we call the accelerated gradient 
method. 

In this paper, we take 



ai 



1 



face 
UiVi 



(59) 



where the factor /acc may be chosen empirically to maxi- 
mize the speed of convergence. We choose /acc = 5. This 
factor seems to take care of the fact that h is smaller 
than h by an order of magnitude. With this choice for 
the scaling factor, aiTii ~ h for deeply bound orbitals 
and aiTLi ~ 5/i for high-lying orbitals. Thus, all the or- 
bitals evolve at roughly the same pace. 
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Incidentally, by using the approximate equality symbol 
in Eq. I|59|) . we mean the utilization of some common 
recipes for numerical calculations like taking a moving 
average and restricting on the maximum values. 

It should be noticed here that one should also mod- 
ify the form of multipliers Ay when one introduces the 
acceleration factors 0;^ 7^ 1. In principle one may reach 
the solution by using Eq. (|49|l because transformations of 
variables do not change the location of the minimum of 
R. Nevertheless, one should modify Eq. (|49|l for a practi- 
cal reason. In calculating the gradient vector whose exact 
form is given by the second member of Eqs. H48|l . the last 
term takes much more computing time than the first two 
terms due to SXjk/Sip*. One can drop the last term if 
the orthogonality relation (|12(l is fulfilled along the path 
of the evolution. Let's suppose that the relation is sat- 
isfied before a gradient-method step is taken and require 
that it is conserved without the dropped term to the first 
order in At after the step, i.e., 

{Am = S,,, (^^IV-;) = <5,, +0 (( AO') , (60) 

with V-'i = ^Ji + Sipi and J^/'i given by Eq. igSJl . Substitut- 
ing Tpi and Tp'^ in Eq. (|60|l and requiring the hermiticity 
(|^ result in 

A., = {ip,\ [aOi, + ajHj) |^/'^). (61) 

This new form of Ay , as well as the original form of 
Eq. (|49|) . fulfills the requirement that it should agree 
with the expression H51|) at the solution. The original 
form differs from the new form, however, before reach- 
ing the solution if ai ^ aj . Consequently, only the new 
form conserves the orthogonality during the course of the 
evolution. 

We have indeed suffered from large errors of orthog- 
onality by using the original form. On the other hand, 
by using the new form, we have observed that the er- 
ror does not grow but decreases without performing ex- 
plicit orthogonalizations. The reason for this stability 
will probably be found in the second order terms in At 
neglected in Eq. (|^ . 

There is an auxiliary method to speedup the conver- 
gence. It is to diagonalize the HFB super matrix given 
by Eq. H38|) in the subspace spanned by iV^ canonical 
orbitals. Then one obtains iVw quasiparticlc states with 
positive energy, whose two-component wave functions de- 
fined in Eq. |SJ are expressed as, 

N„ N„ 

Mr) = J2^jir)U,„ (p,(r)=X]V,(r-)V^,., (62) 

where {Uu, • • • , Um^i, Vu, • • • , VN^i)"^ is the ith normal- 
ized eigenvector of Eq. H38|l . The HFB ground state is 
constructed as the vacuum of these quasiparticles as in 
Eq. Q. By diagonalizing the one-body density matrix 
of this vacuum in this subspace, V*V'^, one obtains a 



renewed set of canonical orbitals as the eigenvectors and 
vf as the eigenvalues. If the initial set of canonical or- 
bitals are taken from the exact solution, the renewed and 
the initial sets are identical. However, if the initial set 
is taken before the state converges to the solution, the 
renewed set corresponds to a better solution than the ini- 
tial set. It is because the above procedure gives the exact 
variational minimum 21, 30] in the subspace spanned by 
the initial canonical orbitals. The gradient-step method 
takes care of both the variation inside this subspace and 
the optimization of the subspace itself. The diagonal- 
ization performs only the former part but perfectly in 
a single step. This diagonalization method is not indis- 
pensable but useful to obtain the solutions. It makes 
the convergence more robust and somewhat quicker if 
it is performed after every ~100 gradient- method steps. 
Its effect seems to saturate with this interval since using 
shorter periods does not lead to noticeable improvements. 
A numerical example is given in Sec. IVI AI 

Incidentally, the idea of this diagonalization method 
originates in the two-basis formalism of HFB Q , in which 
the quasiparticlc Hamiltonian is diagonalized in the sub- 
space spanned by low-lying HF orbitals (eigenstates of 
h), not the canonical orbitals as in this paper. 

The relation between these quasiparticles and the true 
quasiparticles defined in the full space is discussed in 
Sec lyTDl 



VI. RESULTS OF NUMERICAL 
CALCULATIONS 

In this section, we show the results of numerical cal- 
culations using a newly developed computer program 
II Ellil El of the canonical-basis HFB method in a 
three-dimensional Cartesian mesh representation. The 
parameters of the model are as follows: For the mean- 
field interaction, Skyrme SHI force is employed but its 
spin-orbit term and the Coulomb force are not included. 
For the pairing interaction, pc — 0.32 fm^'^ and pc — 00 
are used. The value of Vp is changed depending on fee (= 
1.7 - 5 fm~^) and TV^ (= 14 - 210) to make the average 
pairing gap to be 2 - 2.5 MeV. Unless the values are spec- 
ified, fcc=2 fm-i, Vp = -1050 MeV fm^, and = 21 are 
used. With N„ = 21, the number of orbitals amounts to 
three times as large as the number of the nucleons, since 
each orbital can be occupied by four nucleons due to the 
spin-isospin degeneracy. 

The calculated state is the axially symmetric prolate 
solution of x4Sii4. The single-particle wave functions are 
expressed on a Cartesian mesh having (in the default set 
up) 29 X 29 X 29 mesh points with mesh spacing a=0.8 
fm. As the boundary condition, we assume a cubic box 
whose infinitely high walls are located at the 0th and the 
30th mesh points. The center of mass of the nucleus is 
placed at the center of the box. The 17-point formu- 
lae [sj are employed to calculate the first and second 
derivatives. In applying these formulae, wave functions 
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are assumed to be antisymmetrically extended beyond 
the walls. Spatial integrals are done using the mid-point 
(=trapezoidal) rule. 



A. Convergence to HFB solutions 

In Fig. 12 examples of the convergence curves are shown 
for several versions of the gradient method. The calcula- 
tions are done using the default parameters. The initial 
wave functions of the canonical orbitals are taken from 
the eigenstates of the standard harmonic oscillator poten- 
tial with a quadruple deformation of /9 = 0.3 and 7 = 15°. 
The -D2/1 symmetry is conserved to a very good accuracy 
during the course of the gradient-method evolutions. 

For the size of an imaginary time step, we use fevo = 
1.0, defined by Eq. H5t)|) . for all the curves. This value 
of /ovo roughly optimizes the convergence speed for the 
fastest method. For slower methods, /ovo rnay be slightly 
larger but such a small change does not affect our discus- 
sions below. 

In connection with the step size, we also adopt a 
recipe from the HF-I-BCS method of Ref. accord- 
ing to which the changes of h and h after each gradient- 
method step should be further damped compared with 
the changes of (by a factor of 0.4 in our case). It is a 
precaution against occasional instabilities like oscillating 
behaviors. This additional damping does not affect our 
conclusion on the comparison of the methods, either. 

The top portion of Fig. 13 shows the total energy i^tot 
of Eq. (|22|l . The middle portion shows the maximum 
over i — 1, - ■ ■ , of the error of the second equality of 
Eqs. (|48|l (neglecting the contribution from the error of 
orthogonality). 



(63) 



which is an indicator of the accuracy of HFB solutions. 
As for HE, this quantity is reduced to the maximum over 
« = 1, 



of 



A: 



.HF 



^\h^b^-e^^b,\ = VWh^\W^¥¥\W, (64) 



which is just the energy width of j-i/'i) for h. The bot- 
tom portion shows the mass quadrupole deformation pa- 
rameter f3 defined as the general definition 21] times 
•\/45/167r ~ 0.95 disregarding nucleon form factors. HE 
result for /3, converging to 0.44, is omitted. The abscissae 
are common to all the portions and designate the number 
of gradient-method steps. 

Thin solid lines in all the three portions are the re- 
sult of the nai've (i.e., with ai = 1) gradient method. 
These curves demonstrate that one can indeed obtain an 
HFB solution with the canonical-basis HFB method on 
a mesh, because after a million steps, i?tot and P appear 
to reach a plateau and Ae^^^ becomes as small as 0.1 
keV. We have also obtained similar convergence curves 




FIG 

text 



10-" 10 10" 

GRADIENT STEP 

2: Convergence to HF and HFB solutions for ^®Si. See 
for explanations. 



for other quantities like the axial asymmetry 7 and the 
r.m.s. radius. There is a serious problem, however, of the 
obvious slowness of the convergence. 

Dot curves show the convergence to an HE solution 
with pairing interactions turned off. (With "w/D" in the 
legend in the top portion of the figure we mean that di- 
agonalization of h is done in the subspace of occupied 
jA = 7 HF orbitals to renew the orbitals after every 
100 gradient-method steps. Without this diagonaliza- 
tion, ipi are not HE orbitals and do not satisfy Aef ^ = 0. 
This diagonalization also makes the convergence some- 
what quicker.) One can see that the HE energy reaches 
a plateau by three orders of magnitude faster than the 
nai've HFB method. In the middle portion, while HF can 
achieve precision of 0.1 kcV with only 1400 steps, HFB 
requires 55000 steps for 10 keV, 3 x 10^ steps for 1 keV, 
and 1.5 x 10^ steps for 0.1 keV precisions. 

A method to improve the convergence speed is the di- 
agonalization of the HFB super matrix in the subspace 
spanned by canonical orbitals. Dot-dash curves are 
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obtained by performing this diagonalization after every 
100 steps. One can see that the convergence speed is im- 
proved by an order of magnitude. However, this progress 
is not satisfactory yet compared with the HF curve. 

Dash curves uses the accelerated gradient method (i.e., 
Q!i 1, expressed as "w/A" in the legend), which brings 
about a speedup by two orders of magnitude. A compar- 
ison with dot-dash curves suggests that the variation of 
the canonical-basis subspace itself is a more important 
factor than a minimization inside the subspace. 

Using both the diagonalization and the accelera- 
tion methods leads to thick solid curves (expressed as 
"w/A,D" in the legend). Steep changes in every 100 steps 
are due to the diagonalizations. One can see that these 
curves converge almost as quickly as the dot curves of the 
HF method. These results demonstrate that canonical- 
basis HFB can be solved without very heavy numeri- 
cal computations by adopting the two improvements de- 
scribed above. 



B. Properties of the canonical- basis orbitals 

In this subsection, we use Fig. 13 to discuss the prop- 
erties of the canonical orbitals in connection with the 
mean- field and the pairing Hamiltonians. An analysis 
using these two Hamiltonians leads to a comprehensive 
understanding of the nature of high- as well as low-lying 
canonical orbitals. The calculation has been done for 
^^Si with iVw = 210 canonical orbitals, which is 30 times 
as many as the number of the orbitals below the Fermi 
level. For the combination of this A^w and kc = 2 fm~^, 
we have adjusted the strength of the pairing force to be 
i;p ~ —837 MeV so that the average pairing gap have a 
reasonable size, 2.0 MeV. 

The result shown in Fig. O has been obtained after 
1.2 X 10^ gradient-method steps: Because of the slow 
convergence of very high-lying orbitals (even using the 
accelerated method), we continued the gradient-method 
evolution until no changes are going on in any of e^. The 
last change took place during 5 x 10'*th to 8 x 10'*th steps, 
in which two orbitals gradually increased their energies 
from ~ 50 MeV to 80 MeV. During this change Aef 
increased slightly up to 3 x 10~^ MeV. In the last 4 x lO'* 
steps there were no changes and AeJ^^^ was decreased to 
< 10"^ MeV. 

The top-left panel shows the logarithm of the occu- 
pation probability v"^. The abscissa (common to all the 
four panels in the left column) is the label i of canonical 
orbitals, which are sorted in the descending order of vf. 
First seven orbitals have occupation probabilities close to 
unity. Above them, the decrease of the probability looks 
rather steady. It is in fact a result of a mixing of several 
sequences: The top-right panel plots the same vf data 
versus (only > part is shown), in which one can 
see three sequences. Let us denote the number of the os- 
cillator quanta by A"osc- We have confirmed that the top 
sequence corresponds to A'osc = I subshells (4 < Z < 9 



are completely seen in e > region) . The lower bunch of 
dots are for / < A^osc - 2. For e > 40 MeV, I = A^osc - 2 
sequence is bifurcated from I < A^osc — 4 sequences. These 
observations indicate that the canonical orbitals have the 
same shell structure as that of the harmonic oscillator at 
least up to several tens of MeV. Then, a question arises 
what generates this shell structure. The second and third 
top portions give clues to answer this question. 

The panels in the second top row show the mean-field 
Hamiltonian h and its expectation value and width for 
each canonical orbital. In the right panel, the mass term 
B and the potential term V of h are plotted as functions 
of the distance from the center of the nucleus. Their 
profiles in the symmetry (z) axis as well as those in the 
equatorial (xy) plane are shown. In the left panel, mean- 
field energy of each canonical orbital is designated using 
a dot and a vertical line. The dot corresponds to the 
expectation value while the vertical line connects be- 
tween ei ± Aei with 

e, = {iP,\h\iP,), Ae, = y'(V,,|/i2|^^) _e2_ (55) 

One can see that first seven orbitals, which are below 
the Fermi level, have very narrow widths (30 - 760 keV): 
They are approximate eigenstates of h, i.e., the (occu- 
pied) HF orbitals. Above the Fermi level, the width be- 
comes larger with increasing expectation value and thus 
the wave functions of canonical orbitals begin to diverge 
from those of (unoccupied) HF orbitals. At the crossing 
point from negative to positive Ci, there seems to be only 
a smooth continuation of this trend. This is an essential 
difference between the HFB canonical and the HF or- 
bitals. Concerning the HF orbitals, may be estimated 
in the Thomas-Fermi approximation, where the number 
of levels of h below e is given as, 

The widths of HF orbitals are always zero. The solid 
curve designates the function i = T{e). It is owing to 
the finite volume of the normalization box that r(e) is 
not infinite for positive e. One can see a distinct differ- 
ence between this curve for the HF orbitals and the dots 
for the HFB canonical orbitals: The level density is by 
far sparser in the latter. (It is indeed the reason to use 
canonical orbitals as discussed in Sec.m) It seems diffi- 
cult to relate the level structure of canonical orbitals to 
h. 

Incidentally, one can see a bifurcation of the sequence 
of the dots {i, e^) at z ~ 30. Its origin is the subshell 
structure: Upper sequence corresponds to high-angular 
momentum (I = iVosc) subshells while the lower one to 
I < Nose — 2 subshells. 

The third panels from the top show similar graphs to 
the second top panels but for the pairing Hamiltonian h. 
In the right panel, attention should be paid to the scale 
of the ordinate, which is smaller by an oder of magnitude 
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FIG. 3: Relations between the properties of canonical-basis orbitals and the mean-field and pairing Harniltonians for ^*Si. 
Panels in the left column show occupation probability, mean-field energy, pairing energy, and radial density distribution for 
each canonical-basis orbital. Panels in the right column show the occupation probability (again, but versus the mean-field 
energy), mass and potential terms of the mean-field Hamiltonian, those of pairing Hamiltonian, and the density distribution of 
the nucleus. See text for explanations. 
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than that for h. It is because the expectation value of h 
is nothing but the negative pairing gap (i.e., = ^^i) 
and thus its size is only a few MeV. In the left panel, 
dots are plotted at and vertical lines connect between 
ii ± Aci, where 

{i^,\h\i^,), Ae~ ^ ,J {i^,\h^\^,) - il (67) 

One can see an opposite trend of the width between 
h and h: In the former (latter), the width is smaller 
for lower (higher) expectation values. It suggests a view 
that high-lying canonical orbitals are closely related to h. 
This can be understood in terms of the state dependent 
Hamiltonian, "Hi — vfh+UiVih. Because of this form, it is 
close to h (h) for deeply bound (highly excited) canonical 
orbitals. 

The solid curve designates the Thomas-Fermi estima- 
tion of r(e), the number of levels of h below e. F can be 
obtained by replacing r,e,V,B with T,e,V,B, respec- 
tively, in Eq. (|66|l . One can see that this curve agrees 
quite well with the dots. (Sorting the orbitals more ap- 
propriately in the ascending order of makes the agree- 
ment better.) It supports the view that h roughly de- 
termines the level structure of canonical orbitals. It also 
guarantees that the application of the Thomas- Fermi ap- 
proximation to h provides a rather precise method to 
count the number of canonical orbitals. 

One can see another difference between h and h that 
Ci becomes positive for i > 15 while is still negative 
at i = 210. It means that h has by far more number of 
bound states than h. Its origins in the Thomas-Fermi 
approximation are the larger ratio V/B than V/B inside 
the nucleus and the vanishing behavior of -B as r — > oo. 

It is an interesting question whether there are finite or 
infinite number of canonical orbitals with vf > 0. Non- 
zero occupation probability means > and thus e; < 
0. Therefore one has to count the number of orbitals with 
negative e. To give a Thomas-Fermi estimation of this 
number, let us suppose the form of density tails to be, 

P = Pin , T^fin . (68) 

KT KT 

For the kinetic-energy density inside the nucleus, a rea- 
sonable estimation is fjn — fcpPin, where hk^ is the Fermi 
momentum. Using a relation V^p = K?p and assuming 
p Pc appropriate in the peripheral region, one can 
show that V{r)/B{r) is independent of r. Consequently, 
the number of orbitals becomes proportional to the vol- 
ume of the integral region Vbox- The explicit form of the 
result is 

r(e = 0) -FboxFo, Fo = ^ (2fc2 _ _ fc2^ . (69) 

Therefore, the number of non-zero-occupation canonical 
orbitals is infinite. 

The bottom panels show the radial distributions of the 
nucleon density (right panel) and that of each canonical 



orbital (left panel). The former is shown in two directions 
while the latter is integrated over the angles as, 

p,{r)^r^ J mr)\dr, (70) 

and thus normalized as pi{r)dr = 1. Vertical lines are 
drawn in the interval of r where pi{r) > 0.001 fm^^. The 
width of the lines are proportional to Pi{r), except that 
the width is saturated at 0.8 (in the scale of the abscissa) 
for Pi{r) > 0.5 fm~^. One can see that the canonical 
orbitals are localized in the neighborhood of the nucleus 
and only gradually shifted outward with increasing i. It 
is not due to the finite box size because even the 210th 
orbital is located far from the boundary, which is at 12 to 
12\/3 fm. The location of the maximum of pi{r) agrees 
fairly well with the classical turning point of V for energy 
ii (as seen in the third top panel in the right column). 
This vindicates the view that h determines the spatial 
extent of high-lying canonical orbitals. 

Let us give another fact to confirm this view. Evalu- 
ating Ea. (|69|l using kc — 2 fm~^, k — 1 fm~^ (see next 
subsection), and fcp = 1.4 fm^^, one obtains Fg = 0.14 
fm"^. For r{i = 0) = = 210, the value of Fbox be- 
comes equal to the volume of a sphere of radius 7.1 fm. 
This roughly agrees with the fact that the 210th orbital 
has the maximum density at r = 7.5 fm. 

The discussions of this subsection can be summarized 
into three points: i) Canonical orbitals well below the 
Fermi level are localized by the mean-field potential, ii) 
Highly excited canonical orbitals are localized by the 
pairing potential, iii) The pairing Hamiltonian can have 
infinite number of localized orbitals due to the vanishing 
mass term at large radii. 

It should be noted here that we do not insist that h 
alone determines high- lying orbitals to the details. One 
has to keep in mind an asymmetry between low- and 
high- lying orbitals. Because lower- lying orbitals have 
much stronger influences on the total energy of the nu- 
cleus, the determinations of their wave functions are 
hardly affected by the orthogonality conditions with the 
higher-lying orbitals. However, the opposite is not true: 
High-lying orbitals must be orthogonal to the orbitals be- 
low the Fermi level, which are determined mainly by h. 
Consequently, although the most essential points can be 
explained by the idea of the change of the Hamiltonian 
from h to h, there may still remain something delicate in 
the treatments of high-lying canonical orbitals compared 
with those of low-lying ones. 



C. Tails of canonical- basis wave functions 

Wave functions of bound states have an asymptotic 
form for large radii of 

mr)\-^. (71) 
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FIG. 4: Radial density profiles of canonical orbitals of Si. 
The thick line represents the total nucleon density (p) while 
the ith thin line from the top designates the angle- averaged 
squared wave function of the ith canonical orbital (IV'iP)- The 
canonical orbitals are labeled from 1 to A'w=21 in the descend- 
ing order of v^. The scale in the ordinate applies only to p 
and IV'iP- The other thin lines are shifted downward so that 
they do not overlap with each other. 



In HF, negative energy orbitals have the form (|71|) with 
Ki determined by as fiKi = \J —Imti. In HFB, the hole 
component of quasiparticle wave functions (}f>) has this 
form with k determined by the excitation energy E'qp as 

As for the canonical orbitals of HFB, however, one can- 
not relate Ki to any sort of energies in a simple manner 
due to the complex influences of the requirement of or- 
thogonality, which is rooted in the lack of a common 
Hamiltonian to all the orbitals. The only possible state- 
ment is that it is not larger than the smallest value of 
K for quasiparticles, whose lower bound is given by the 
Fermi level as in Eq. (|14|l . 

In this subsection, we compute Ki of canonical orbitals 
numerically using angle-averaged density: |?/'i(r)p — 

J IV'i(''')P'^^ ■ III order to see the behaviors at large 
radii, we have used a large box containing 57 x 57 x 57 
mesh points. The edge of the box is 58a=46.4 fm. The 
nearest walls from the center of the nucleus is at a dis- 
tance of 23.2 fm and the farthest (i.e., the corners) at 
40.2 fm. For the sake of using the large box, we restrict 
the number of orbitals to iVw = 21, for which an ade- 
quate strength of the pairing interaction is Vp = —1050 
MeV. 

Figure ^ shows |i/'i(r)p of canonical orbitals, all of 
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FIG. 5: Coefficient Ki for the tails of canonical-orbital wave 
functions versus the expectation values of the mean-field 
Hamiltonian e^. The orbitals are the same as those shown 
in Fig. 4. Open circles and crosses are for k at r = 6 fm and 
r = 20 fm, respectively. 

which have exponentially damping tails. This result 
demonstrates that the canonical-basis HFB method can 
convert Gaussian tails of initial wave functions, for which 
we use eigenstates of a harmonic oscillator, into the cor- 
rect exponential form. 

It is evident that the first orbital changes its slope at 
around r = 10.5 fm. By looking at the other curves 
minutely, one can also find such changes in several other 
orbitals. Therefore, we give two discussions, first on the 
behaviors near the nucleus and second on those far away. 

Figure shows Ki as a function of e^. Open circles 
denote Ki determined by the logarithmic derivative of 
r|i/'i(?')| near the nucleus (r = 6 fm). One can see 
less than 21 circles because several pairs of circles over- 
lap exactly due to the degeneracy for the sign of Iz {z- 
component of the orbital angular momentum) . The dash 
and dot curves represent functions JiKi = \J ~2mti^ and 
fiKi = \J'2,m{ti — 2eF), respectively. The open circles 
seem to agree fairly well with a function 

hK, = V2m(|e, -erl -ep), (72) 

which coincides with the dash curve for < ep and the 
dot curve for ti > ep. 

By noting that the dash curve expresses k for bound 
HF orbitals, one may understand the < ep part of 
Eq. H72() as indicating that canonical orbitals below the 
Fermi level have approximately the same tails as HF or- 
bitals. 

An alternative interpretation is also possible which 
refers to the full domain of Eq. (|72l) : Canonical orbitals 
with Ei — epdzEqp have the same k as the hole component 
(/? of a quasiparticle state with Eqp > 0. It may suggest 
that canonical orbitals below (above) the Fermi level are 
related to the hole components of hole-like (particle-like) 
quasiparticle excitations. However, in order to confirm 
it, we need quasiparticle states in the full space, which 
are outside the scope of this paper. 

Cross symbols plotted in Fig. |5l are Ki calculated at a 
farther radius (r — 20 fm). In this case, concerning the 
deepest and highest orbitals, k drops toward the min- 
imum value at = ep. This may be ascribed to the 
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fact that the canonical orbitals are constructed by mixing 
up the hole part of all the quasiparticle states and thus 
the component having smallest k dominates in all the 
canonical orbitals at large radii. Probably, all the canon- 
ical orbitals have the same k (— y^2m(Eqp^min — ^f)/^ 
where -Eqp^min is the smallest quasiparticle energy) at suf- 
ficiently large radii. 



D. Quasiparticle states 

Let us discuss on the quasiparticle states defined by 
Eq. in the subspace spanned by iV^ canonical or- 
bitals. These quasiparticles and the true quasiparticles 
defined in the full space are quite different. The former 
are only a set of tools for the minimization in a small 
subspace, while the latter have physical information on 
the excitation modes of the nucleus. 

The width AEqp (in the full space) of the excitation 
energy Eqp of the ith quasiparticle (in the subspace) is 
given by 

(Ai?W)' = I (|0Ur)P + |^Kr)P)dV-(£;«)', 

(73) 

\ h -h + eFj\ip^J' 

where the operation of h and h should be evaluated in 
the full space (i.e., in the mesh representation), not in 
the A^w-dimensional subspace. 

Quasiparticles can be characterized by the norms of the 
particle and hole components as well as by the excitation 
energy. These norms are defined as 

Uf^J \Mr)fd\, = J \^^ir)fd\, (74) 

for the particle and hole components, respectively. They 
are normalized as C/f + Vj^ = 1 according to Eq. jnj. 
Quasiparticle states are called to be particle-like (hole- 
like) ifUfiV,')>^. 

FigureElshows Eqp and Ai?qp for A^w = 210 (with Vp = 
—837 MeV fm^). One can see that hole-like quasiparticles 
have very small widths. Even the deepest hole state (Is) 
has a relatively small width (1.7 MeV). On the other 
hand, the widths of particle-like quasiparticles amount 
to a few tens of MeV: They are very different from the 
true quasiparticles defined in the full space. 

For the sake of comparison, quasiparticle levels ob- 
tained in a smaller subspace {N^ — 21, Vp — —1050 
MeV fm^) are also shown in the top-right window. By 
comparing the results of two calculations, one can see 
that increasing the number of canonical orbitals leads 
to the addition of high-excitation particle-like quasipar- 
ticles, while it does not change very much the hole-like 
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FIG. 6; Excitation energies Eqp (dots) and their widths 
A£'qp(half of the length of vertical bars) of quasiparticle states 
of ^*Si plotted versus the norm of the hole component. Cal- 
culations are done using the HFB solution for A'w=210 shown 
in Fig. 3. The top-right window shows the result for N^=21 
for comparison. 

states and low-excitation particle-like states. Thus, en- 
largement of the canonical-basis subspace is a very inef- 
ficient way to decrease the energy widths of particle-like 
quasiparticle states. 

The differences from true quasiparticles become more 
evident by looking at the wave functions. For Eqp > — A, 
the particle component (f) should be an oscillating func- 
tion of r 3]. However, (p of quasiparticles defined in the 
canonical-orbital subspace has an exponentially damp- 
ing tail. It is only natural because they are expressed as 
linear combinations of such functions. 

It should be stressed again that the quasiparticles in 
the canonical-orbital subspace are nevertheless a useful 
auxiliary tool to obtain the HFB ground state. How- 
ever, if one needs true quasiparticle excitation modes, 
the quasiparticles in the subspace are useless. One has 
to calculate with some other method the true eigenstates 
of the quasiparticle Hamiltonian. A good news is that the 
Hamiltonian (being composed of h and h) is already ob- 
tained and thus iterations for self-consistency are not nec- 
essary. This will make the numerical calculations rather 
inexpensive. 



E. Choice of the value of kc 

Our canonical-basis HFB method has several parame- 
ters to specify the pairing force, i.e., the over-all strength 
Vp, the parameter for the momentum dependence 
the number of canonical orbitals having non-zero occu- 
pation probabilities A^w, and parameters for density de- 
pendences. However, available experimental data do not 
provide sufficient information on the pairing correlation 
for precise determinations of so many parameters. There- 
fore we tune only Vp precisely so that the resulting pairing 
gap has a desired value, while assuming some reasonable 
values for the other parameters. However, this does not 
guarantee that their values may be chosen completely ar- 
bitrarily. In this and next subsections we discuss on the 
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FIG. 7: Dependence of the average pairing gap A on the 
number of canonical orbitals N„ for three sets of pairing in- 
teraction parameters kc (fm~^) and Vp (MeV fm'^). See text 
for explanations. 



choices of the values of kc and N^,, respectively. 

Pairing-channel interactions were taken into account 
when Skyrme forces SGII 33| and SkP 0] were deter- 
mined. Their values of fee {— •\/— to(l — xo)/ti{l — xi)) 
are 2.5 fm~^ and 4.3 fm"^, respectively. From this result, 
values from 2 to 5 fm~^ seem equally reasonable. 

From a pragmatic point of view to use the force in the 
canonical-basis method, kc must be smaller than 6 fm~^ 
in order to avoid the point collapse (see Appendix IB|| . 
Values smaller than 1.5 fm~^ are not preferable, either, 
because such small kc makes the solutions very difficult 
to obtain. 

If Vp is adjusted for each value of kc to reproduce the 
same value of the average pairing gap, which we define 
El as , 



N„ 



i=r 



1=1 



(75) 



resulting properties of the nucleus are almost indepen- 
dent of kc except the pairing density, which becomes more 
diffuse when smaller kc is used. Unfortunately, pairing 
density is rather difficult to determine experimentally. 

Incidentally, the properties of each canonical orbital 
are rather sensitive to kc'. and {rf) increase while 
decreases with decreasing kc. Only stays roughly con- 
stant. Nevertheless, properties of the nucleus (except the 
pairing density profile) are hardly changed, as far as the 
average pairing gap is kept unchanged. 

From the allowed interval 1.5 fm~^ < fee < 6 fm^^, we 
have chosen 2 fm~^ in this paper. There is no physical 
reason for this choice. However, smaller values of kc lead 
to a favorable situation that the dependence on be- 
comes weaker and the determination of A^w can be less 
precise. Fig. shows the dependence of A of Eq. iff^ 
on iVw. Sohd circles, open circles, and triangles are ob- 
tained with kc= 5, 3, and 1.7 fm"-'^, respectively. The 
values of are given in the legends of the figure in units 
of MeV fm'^. Here, Dp is not changed as a function of 
but is determined for each kc such that the resulting A is 
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FIG. 8: Pairing force strength Vp adjusted to give A = 2.5 
MeV as a function of the number of canonical orbitals A^w for 
^®Si. Dots denote the strengths obtained by the adjustment 
while the curve represents a function fitted to the dots. 
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FIG. 9: Dependence of the total energy Etot on the number 
of canonical orbitals A'^w for '^^Si. See text for explanations. 



roughly the same at ^ 28. One can see in the figure 
that, for kc = 1.7 fm"-'^, dA/dN„ is almost zero around 
~ 80. The emergence of this plateau is explained by 
the vanishing of the pairing interaction at relative mo- 
mentum kc ■ For much larger N^, A will increase again 
infinitely 



F. Dependences on Aw 

We next study the dependence on the number of canon- 
ical orbitals N^. Here we adjust the strength Vp for each 
value of A^w to keep A constant. The used values of Vp 
are shown in Fig. |S| 

In Fig. 1^ we show the total energy of ^^Si as a func- 
tion of A^w. It is obtained as follows: Starting from har- 
monic oscillator eigenstates, we obtain an HFB solution 
for — 70. Then, we remove the canonical orbital 
with the smallest v'f and again obtain the HFB solution 
for A'w = 69 by the gradient method. By repeating this 
procedure, we obtain solutions for > 14. This figure 
shows that the magnitude of the influence of A^w on the 
total energy is 300 keV or less when Vp is determined for 
the average pairing gap. Other bulk properties of the nu- 
cleus are also roughly independent of A^^: Deformation 
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FIG. 11: Dependence of the average pairing gap A on the 
mass number A. The strength of the pairing interaction is 



fixed at 



4000 MeV fm"' while is changed in two 



ways as shown in the legend. 



culations, the normalization box is gradually expanded 
with increasing A: The edge of the box is at least 2.5 
times as large as the liquid-drop-model diameter. We fix 
Vp at —1000 MeV fm'^. Instead, we change versus A. 

When we assume a relation AN^, = 2A, we cannot re- 
produce the empirical trend that A is a decreasing func- 
tion of A. We can improve the result by using different 
relations, e.g., ANy^ — A + 125. However, the goal, a uni- 
versal pairing force, should be more than an empirical 
formula expressing N^j as a function of A. It will be one 
of our future challenges. 



FIG. 10: Dependence of the mean-field energies of the 
canonical orbitals on the number of the orbitals A^„ for ^®Si. 
See text for explanations. 



parameters /? and 7 are within 0.374 ± 0.01 and ±0.3°, 
respectively. R.m.s mass radius is 3.201 ± 0.004 fm. 

Such insensitivity can also be seen in the properties 
of individual canonical orbitals. In Fig. ^1 mean-field 
energies of canonical orbitals are plotted versus N^,. 
One can see that below 5 MeV stays almost constant. 
Higher levels are not so constant but only slightly de- 
creasing. We have also confirmed that other quantities 
also remain roughly constant. 

Therefore, unlike kc, the effects of changing Nyj is al- 
most canceled by the adjustment of Vp. 

G. Application to heavy nuclei 

Finally we demonstrate the applicability of the 
canonical-basis HFB method to nuclei heavier than ^^Si. 
We have encountered no special difficulties at least up to 
A = 252 except the increase in the computation time per 
gradient-method step. 

For example. Fig. II II shows the dependence of the av- 
erage pairing gap A on the mass number A. In the cal- 



VII. CONCLUSIONS 

In this paper we have presented a method to obtain 
solutions of the HFB equation expressed in the canonical 
form, i.e. in terms of BCS-type wave functions, without 
using quasiparticle states. The method is fit to three- 
dimensional coordinate space representations and is ad- 
vantageous to describe simultaneously arbitrary defor- 
mations, long and short density tails, and pairing cor- 
relations involving states in the continuum of the HF 
Hamiltonian. 

We have improved the speed of the convergence to HFB 
solutions by a few orders of magnitude. For this purpose, 
we have modified the gradient method and derived the 
appropriate form of the Lagrange multiplier functionals 
for the orthogonality between canonical orbitals. The 
two-basis method has also been adapted by replacing the 
HF basis with the canonical basis to further speedup the 
convergence. 

For delta-function pairing interactions, the wave func- 
tions of canonical orbitals have been shown to shrink in- 
finitely into a spatial point once their occupation proba- 
bilities become less than some critical value. This prob- 
lem is peculiar to the canonical-basis HFB method. To 
avoid it, a repulsive momentum dependent term has been 
added to the pairing interaction. 
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As a byproduct of the introduction of this term, the 
nature of high-lying canonical orbitals has been greatly 
clarified as approximate eigenstates of the pairing Hamil- 
tonian. We have demonstrated that the level density and 
the spatial extent of high-lying canonical orbitals can be 
computed by applying the Thomas-Fermi approximation 
to the pairing Hamiltonian. 

The obtained canonical orbitals have been examined 
in many aspects. One of the results is that canonical or- 
bitals closer to the Fermi level have wave functions with 
longer tails. At very large radius, however, all the or- 
bitals seem to change the slope of their tails to that of 
the orbital closest to the Fermi level. Another result is 
that quasiparticle states cannot be efficiently expanded 
in the canonical basis. Therefore the canonical orbitals 
perfectly describing the ground state are not sufficient to 
treat excitations. 

The pairing force strength has been adjusted to repro- 
duce the average pairing gap as a function of the number 
of canonical orbitals, which is used instead of an energy 
cut-off. Properties of the nucleus have been found to be 
almost independent of the number of canonical orbitals 
and rather insensitive to the relative strength of the mo- 
mentum dependent to independent terms of the pairing 
force. 

Although most of our discussions refer to calculations 
for ^^Si, we have also checked the applicability of the 
method to heavy nuclei up to A— 252. Checks have also 
been done concerning the precision of the mesh represen- 
tation in connection with this method. 

We believe that the canonical-basis HFB method has 
the potential to become the standard approach to treat 
neutron-rich nuclei in HFB. As a next step, we are going 
to extend the formalism and the computer program con- 
cerning the treatment of N ^ Z nuclei and the inclusion 
of the spin-orbit and the Coulomb interactions. 



APPENDIX A: PRECISION OF THE MESH 
REPRESENTATION 

In this Appendix, we examine the precision of the 
Cartesian mesh representation when it is applied to the 
canonical-basis HFB method. Discussions are given con- 
cerning the dependencies on i) the mesh spacing, ii) the 
approximation formulae for the derivatives, and iii) the 
size of the normalization box. The pairing force strength 
used is Up = -1032 MeV fm^ with N„ = 21. The state 
is the prolate ground state of ^*Si. 

i) We discuss the dependence on the mesh spacing a. 
For this purpose, we calculate the HF and HFB solu- 
tions for various value of a while keeping the size of 
the box constant. For the discrete approximation of the 
first and second derivatives, we use the 17-point formula 
[Slf. According to the calculation using the finest mesh 
(a = 0.3 fm), Etot = -257.723 MeV, A = 2.424 MeV, 
and P = 0.3766. 

The mesh is changed in the following manner: We 
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FIG. 12: Errors of the total energy Etot (top portion) and 
the quadrupole deformation parameter /3 (bottom portion) 
due to the finite mesh spacing a for the prolate solution of 
^*Si. The result of HF (HFB) calculations are designated 
by triangles (circles) connected with dot (solid) lines. In the 
bottom portion, triangles and circles are open (solid) when 
the error is positive (negative). 



use TVmcsh mesh points in x, y, and z-directions. Then 
the length of the edge of the normalization box is L = 
(TVmesh + l)a. We fix L at 24 fm so that increasing A^mcsh 
from 15 to 79 corresponds to decreasing a from 1.5 fm 
to 0.3 fm. The error is defined as the deviation from the 
result with a=0.3 fm. The center of mass of the nucleus 
is placed at the center of the box. 

The top portion of Fig. [T^ shows the error of the total 
energy i^tot as a function of the mesh spacing a. The er- 
ror of the corresponding HF calculation is also shown 
for comparison. Before obtaining this result, we sus- 
pected that HFB solutions suffer from larger errors than 
HF solutions because the former involve states above the 
Fermi level while the latter depend on only levels below 
it. Higher levels contain more amount of high momen- 
tum components which are less accurately treated with 
a finite mesh spacing. However, this figure shows that 
the errors are of the same order. The principal reason 
may be that errors from high-lying levels are reduced by 
factors vf or UiVi. 

The bottom portion of Fig. 1121 shows the error of the 
quadrupole deformation parameter (3. Unlike i?totj the 
error of P takes both positive and negative values. For 
this quantity, the HFB errors are several times as large 
as the HF errors on the average. This may be related to 
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FIG. 13: Dependence of the error of the HFB total energy 
Etot of ^*Si on the discrete approximation formulae used to 
evaluate the derivatives. Three curves correspond to different 
values of the mesh spacing a. 
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FIG. 14: Errors of the total energy Etot (on the left-hand 
side) and the r.m.s. mass radius rrms (on the right-hand side) 
due to the finite box size L. 



the strong influence of the pairing correlation on defor- 
mations. 

ii) We compare some approximation formulae given in 
Ref. lU for the evaluation of the derivatives. Note that 
the Fourier transformation method (or the Lagrange- 
mesh method 33] applied to uniform-spacing meshes) as- 
sumes a different boundary condition from the other for- 
mulae. It assumes the periodic boundary condition and 
the edge of the box is also slightly shortened to A^mcsha- 
The error is defined as the deviation from the most pre- 
cise result, which is obtained with the 17-point formula 
and 0=0.3 fm. We have not chosen the Fourier method 
for this purpose because of a large error due to numerical 
cancellations for more than several tens of mesh points. 

In Fig. 1131 the error of the total energy is shown for 
each approximation formula and for three values of the 
mesh spacing a. This figure tells, e.g., if one uses the 
17-point formula with a — 0.8 fm for ^^Si, the error of 
the total energy is 70 keV. 

One can observe two trends: First, for all the three 
mesh spacings, the error decreases as the formula in- 
volves more number of mesh points (Fourier transforma- 
tion method uses all the mesh points in the box). Sec- 
ond, for every formula, the error decreases when the mesh 
spacing is decreased. Therefore, one can decrease the er- 
ror of the total energy either by employing a more precise 
formula or by diminishing the mesh spacing. 

On the other hand, the other quantities cannot always 
be improved significantly by using more precise formu- 
lae. For example, with a=1.0 fm, the error of the pairing 
gap A does not continue to decrease but seems to almost 
saturate around 10 keV for more than 9-point formu- 
lae (including the Fourier transformation). We have also 
found a similar saturation phenomenon in the error of 
the deformation parameter /3, which is affected strongly 
by the pairing gap. The reason of this saturation is prob- 
ably that the pairing gap is more sensitive than the total 
energy to high-lying states having large momentum and 



thus the lack of momentum larger than n/a in the mesh 
space becomes the main source of the error, while the 
change of the formulae can only improve the treatment 
of momentum components smaller that ir/a. 

iii) We examine the errors due to the finite size of the 
normalization box. We consider ^*Si, which has e-p = 
— 11.6 MeV. Naturally the errors depend on the Fermi 
level. In this paper, however, we can treat only N ^ Z 
nuclei, all of which have roughly the same Fermi level. 

We control the size of the box L in terms of A^mosh while 
fixing a at 0.8 fm. We define the errors as the deviation 
from the result obtained with large L (~ 60 fm). Here, 
the comparison should be made within the same parity 

ofiVmesh 

On the left- and right-hand sides of Fig. ^1 we show 
the errors of the total energy Etot and the r.m.s. radius, 
respectively. We have used L = 24 fm in most of the 
calculations shown in this paper. These figures confirm 
that this value of L leads to sufficiently precise results for 
both quantities. By comparing the errors of HF (crosses) 
and HFB (circles) solutions, one can see that the effects of 
the finite box size are of the same order between the two 
methods. This is because the density tail is determined 
by the Fermi level, which is roughly the same whether 
there exist pairing correlations or not. 



APPENDIX B: THE MECHANISM OF THE 
COLLAPSE OF A HIGH-LYING CANONICAL 
ORBITAL INTO A SPATIAL POINT AND ITS 
REMEDIES 

In this appendix we study the mechanism of "point 
collapse" , a phenomenon that the wave function of a 
canonical orbital shrinks infinitely into a point in the co- 
ordinate space once its occupation probability becomes 
less than some critical value. It is a problem peculiar to 
the canonical-basis HFB method. For this purpose, we 
consider only one canonical orbital explicitly while repre- 
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senting all the others in terms of uniform densities, po, po, 
To, and Tq. The explicitly considered orbital is assumed 
to have a wave function whose spatial part is a Gaussian 
wave packet parameterized with a size parameter w as, 
/ 9 \ 3/4 

^^{r)^U\ w'^/'e-^^/-^\ (Bl) 



The r.m.s. values of the coordinates and the wave num- 



bers of this wave packet are, respectively, y/ (x^) = y/Jy^ 

1 



= ^ = f and /(fcl) = y'(fc2) = y(fcf) = 
can express the densities of the total system as, 

P = Po + Pi, " " ""'^ 

T = To + Tl , 

P = Po + P~l, 

^-tq + ti, fi = 4to |V7/'^(r)|^ , 

where v'^ is the occupation probability of this orbital and 

u = \[\-~v^ . 

The change in the total energy due to the presence of 
the explicitly considered orbital is given by 

[H(r)-Ho(r)]r2dr, (B6) 



Pi =4z;2|V'„(r)r 
Tl =4^;2|VV'^(r)|^ 

Tl 



One 

(B2) 
(B3) 
(B4) 
(B5) 



AStot = 47r 



where H, has been defined by Eq. I|23|) . while Tio by re- 
placing p, T, p, and f in the equation with po, tq, po, and 
To, respectively. For fcc = Pc = oo, one can obtain the 
following expression for Ai^tot (a = 1 is assumed for the 
density dependence of the Skyrme force): 
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FIG. 15: Change in the total energy AiStot due to the ad- 
dition of a canonical orbital with occupation probability 
whose wave function is a Gaussian wave packet parameter- 
ized with a width parameter w. Ai5tot is plotted versus w for 
several values of . Local maxima and minima of the curves 
are designated with circles and crosses, respectively. 



Among the nine coefficients, Ci, C2, C4, and Cg are nega- 
tive (at low density places) while the others are positive. 

Fig. El shows Ai?tot as a function of w for several 
values of u^. In the calculations, we employ a pairing 
force having Vp = —450 MeV fm^^, pc = 0.32 fm~'^, and 
Pc = kc = 00. This strength Vp would give A = 2 MeV 
for A ^ 28 system with = 21 if point collapses would 
not occur. The effects of the momentum and the pairing 
density dependent forces are discussed later. The mean- 
field force is SIII without the spin-orbit term. A = 00 
is used in Eq. (|21|l for the nucleon reduced mass. For 
the background densities, we employ po = Po = lxlO^'^ 
fm~'^ and to = f = 3.6 x 10~^ fm~^, which are typical 
values at peripheral regions where point collapses occur 
in actual HFB calculations for finite nuclei. (A Fermi 
gas relation t — |(-^)2/'^p^/-^ is assumed to determine 
T from p.) 

One can see from Fig.^jthat the total energy becomes 
minimum at the bottom of a dip for small values of u^. 
The dip emerges at w = 0.7 fm when u^— g.7 x 10^^ and, 
as v"^ is decreased, its depth increases while the location 
of the bottom moves to smaller values of w. 

This minimum in the dip corresponds to HFB "solu- 
tions" in finite nuclei reported in Refs. Q| , in which 
the wave function of a high-lying orbital has shrunken 
into a mesh point. HFB states having physically reason- 
able wave functions turn themselves into such strange 
states in the following abrupt way: In the course of 
gradient-method evolutions, suddenly after an almost 
stationary situation has continued, e of a high-lying or- 
bital soars up to a few hundred MeV and simultaneously 
the total energy decreases by a few MeV. 

The dip is caused by the zero-rangeness of the pairing 
force and thus does not correspond to physically mean- 
ingful states. Physical solutions should be found at w > 6 
fm, where A_Etot is slightly negative. Therefore one needs 
to eliminate the effects of such dips in order to obtain 
well-behaving HFB solutions in the canonical-basis rep- 
resentation. 
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It is worth pointing out here that a wave packet with 
w=0.3 fm can be expressed fairly precisely on a mesh 
with spacing a=0.8 frn. It is because the maximum wave 
number expressible with the mesh i s ^ a nd thus the mini- 
mum size of the wave packet of Eq. IjBip can be estimated 
roughly by a relation = :^ ~ f i which leads to 

w - = 0.26 fm for a=0.8 fm. Therefore collapses 
to wave packets as small as w = 0.3 fm are a practical 
problem even when a~0.8 fm. 

Before discussing how to remove the dip, let us clarify 
in an analytical way the mechanism which gives rise to 



the dip. We consider the limit of w 



0, in which the 



nature of the dip becomes most manifest. Therefore, the 
following considerations are for the untruncated Hilbert 
space, not for mesh representations. 

In order to study this limit, we assume a scaling be- 
havior V — jiw^ where w is the location of the bottom of 
the dip. /i and u are positive constants to be determined. 
One can easily confirm that, for v > |, all the nine terms 
in the right-hand side of Eq. ljB7|) have positive powers 
of w and thus converge to zero as w ^ 0. Consequently, 
it must be satisfied that v < ^ for a point-collapse solu- 
tion to exist. For < < |, Ai^tot = C's^i'' {'t:)^ '^'' + 
(lower order terms in and thus the energy diverges 
to -|-oo as w ^ 0. The only possibility for a point- 
collapse solution to exist is in the case oiv = |, for which 
AE^tot = Cifx^ + Csfi'^ + (terms with positive powers in 
w). From a condition dAEtot/d^ = 0, one immediately 
finds that AE'tot becomes minimum at fj? ~ —Ci/2C%. 
The minimum value in the limit w — s- becomes 



4C8 



64 



-VpPc 



1 - 



(B8) 



Its numerical value is —11.7 MeV for the parameters 
used. The essential points clarified here are i) the mini- 
mum of the total energy is realized in the limit w — )■ 
with V oc ufl/'^ and ii) the limiting value of the minimum 
energy is finite. 

One can also derive analytically the location w of the 
barrier top which divides the physical and unphysical re- 
gions. From the condition. 



dw 



w 



0{v^)^Q, (B9) 



one can simply obtain the location of the local maximum 
as, 



3C4 
2C^ 



(BIO) 



by neglecting terms of 0{v^). Note that it is independent 
of V (if V is sufficiently small). With the parameters used, 
w = 1.0 fm. Indeed in Fig. ^1 the location of the top 
of the barrier (designated with open circles) seems to 
converge to ~ 1 fm as is decreased. The height of the 
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FIG. 16: Effects of changing kc on the total energy AiStot 
(left-hand portion) and on the width parameter w of the 
added orbital at local minimum and maximum of Ai?tot. The 
dot, triple-dot-dash, dot-dash, dash, and solid curves are ob- 
tained with kc = 15, 30, 100, 500, and 00 (fm~^), respectively. 
The other parameters of the model are the same as those used 
for Fig. 15. 



maximum is Civ + 0{v'^), which is 0{v) and converges to 
zero as w ^ 0. Therefore, for sufficiently small values of 
V, HFB solutions easily fall into the unphysical dip since 
the mesh spacing is much smaller than I.Ott fm. 

Now, we examine two types of pairing interactions 
which can remove the unphysical dip. 

The first type of the interaction is the momentum de- 
pendent term of Eq. (|19|l . With kc < 00, but with 
Pc = 00, Aii^tot becomes to have additional terms as. 



A£:tot + Ciouv -f Cii^ -I- C12- 



(Bll) 



where Ai^tot is given by Eq. (|B7|) and 



Cio — 
Cii - 
C12 — 



V- 

2k, 



P o ~ 

30 



Since Vp < 0, Cio, Cn, and C12 are positive. The dom- 
inant term in the limit w — > is the one having coeffi- 
cient Cii (C12) for V > (<)3, which behaves as fiw"^^ 
{lJ?"uP'^~^). These two terms overwhelm the negative 
terms in Ai?tot- Consequently wave functions cannot 
shrink to a point. 

However, a dip can still emerge with finite v. Fig- 
ure El shows how the local minimum (the bottom of the 
dip) and the local maximum (the top of the barrier) of 
AE^tot change versus . The left-hand portion shows 
the depth and the height (Ai?tot), while the right-hand 
portion shows the locations {w). 

The solid curves are obtained with kc = 00. The upper 
solid curve stands for the height of the local maximum, 
which converges to zero as v'^ 0, as already discussed 
analytically. (This curve is almost unchanged when kc 
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is decreased to finite values.) The lower solid curve cor- 
responds to the depth of the local minimum, which con- 
verges to —11.7 MeV, in accord with Eq. 1)B8|I . 

With finite values of fcc, the minimum becomes shal- 
lower as is decreased and it disappears at some small 
value of . This disappearance is clearly seen in the dot 
curve (corresponding to fcc=15 fm~^) in the right-hand 
portion of the figure. A numerical search has shown that 
the dip exists only for fcc > 10 fm~^. 

The result with this simplified model roughly agrees 
with realistic HFB calculations: We have observed that 
point collapses can occur for fee > 6 fm~^ in calcula- 
tions like those shown in Fig. [T] The difference between 
6 fm^^ and 10 fm^^ can be partially attributed to the 
finite-point approximation to the derivatives, which un- 
derestimates the kinetic energy by 30% at fc = -, as 
shown in Fig. 20 of Ref. 

Incidentally, the requirement of orthogonality among 
the canonical orbitals does not seem to change the con- 
clusion drawn with this model. The first reason is that 
point collapses in light nuclei are not affected by orthogo- 
nality if the number of canonical orbitals iV^ is relatively 
small. This is because of the symmetry (e.g., Dih) of 
the solution: If a high-lying orbital is the only one in 
an irreducible representation of the symmetry (e.g., in 
Px—— sector when a;-axis is the shortest axis), it is auto- 
matically orthogonal to the other orbitals. In this case, 
it shrinks not to a point but to many points (e.g., two 
points at X = ±-R, y = z = 0, i? '--^ nuclear radius, 
with opposite-sign amplitudes). Indeed, we have experi- 
enced that point collapses occur only when is small 
for relatively small values of k^- The second reason is 
an empirical fact that, with large values of fcc, many or- 
bitals collapse one after another or simultaneously even 
for large N^: The process of point collapses under or- 
thogonality condition may not be so complicated as it 
sounds in such cases. 

The second type of the interaction we examine is the 
pairing density dependent term of Eq. H19|l . With < 
oo, but with fcc = oo, Ai?tot has additional terms to 
Eq. jBTl) as. 



where 



AEtot + Ci3uv + Ci4- 
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All the added terms are positive. Assuming v = pw'^ , 
terms with coefficients C15 and Cig have negative powers 
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FIG. 17: Same as in Fig. 16, but for the effects of pc. The dot, 
triple-dot-dash, dot-dash, dash, and solid lines are obtained 
with pc = 1.5,6,25,100, and 00 (fm~'^), respectively. The 
other parameters of the model are the same as those used for 
Fig. 15. 



of w at = |. Consequently point collapses can be 
avoided by using finite values of pc- 

Figure [T7I shows extrema of AE^tot versus v'^. The solid 
curve corresponds to pc = c» and is the same as that in 
Fig. El The other types of curves are for finite values 
of Pc. One can see that AEtot at the local minimum 
is already very small with pc = 1.5 fm~"^ (dot curve), 
which is much larger than the typical values of p and 
thus leads to only physically acceptable small amount 
of modification of the pairing interaction. This figure is 
consistent with the results of actual HFB calculations: 
From Fig. 4 of Ref. IT3|, one can see that point collapses 
can be avoided with p^ < 1.0 fm~^. 

Mathematically speaking, however, the local minimum 
continues to exist for arbitrary v, which is a different be- 
havior from that with finite fcc- In the limit of w — > 0, 
this minimum behaves as w and AE['^^ — 0. 
Because the depth converges to zero from below, the 
point-collapse limit is unstable and is not very mean- 
ingful. So, let us give only briefly the result of an ana- 
lytical consideration. Assuming v — pw'^ , one can show 
that the local minimum continues to exist in the limit 
w with = 3 and p determined by a cubic equa- 
tion Ci + Ci3 + 2C4A< + 4Ci6M^ = if Ci Cia < 
(i.e., Pc > y/2pQ for po ^ pc). At the local minimum, 
AE['^^ (X V and w oc w^/'^. One can see this scaling be- 
havior with = 3 on the right-hand side of Fig. [T7I 

In the present paper we have used the momentum de- 
pendent interaction because it has a physical origin, i.e., 
the finite-rangeness. However, the pairing density de- 
pendent interaction may also be useful under different 
circumstances because the momentum dependent inter- 
action leads to a divergence of the pairing energy in un- 
restricted Hilbert space [S^: The assumption of > 
made in this paper is not adequate for orbitals with 
k ^ kc . 

To summarize, assuming a wave function of Gaus- 
sian form, we have clarified the mechanism of the point- 
collapse phenomenon occurring in the canonical-basis 
HFB method due to delta interactions. We have also 
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demonstrated the successes of momentum or pairing den- sity dependent interactions in avoiding the collapse. 
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